##### Models with CCPI sample #### 

library(tidyverse)
library(broom)
library(fixest)
library(ggthemes)

source("code/01-importdata.R")
source("code/02-mainresults.R")

ccpi_sample <- read_csv("data/ccpi_dummy.csv") %>% 
  select(iso3c, ccpi_dummy)

analysis_ccpi <- full_join(analysis_df, ccpi_sample, by = "iso3c") %>% 
  filter(year>=2008&year<=2017,
         ccpi_dummy == 1)

## Models with country-specific time trends

### Annual count of disasters, for three lags
m_ccpi_annual1 <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                          scaled_elec_renew, scaled_fit, scaled_quota,
                          scaled_delegates, scaled_gcg,
                          scaled_cf_total_provided, scaled_cf_principal_provided,
                          scaled_cf_total_received, scaled_cf_principal_received)
                        ~ annual_average_temp_lag1 + annual_W_lag1 |
                          iso3c[year] + year, # trends
                        data = analysis_ccpi) %>% 
  as.list()

m_ccpi_annual2 <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                          scaled_elec_renew, scaled_fit, scaled_quota,
                          scaled_delegates, scaled_gcg,
                          scaled_cf_total_provided, scaled_cf_principal_provided,
                          scaled_cf_total_received, scaled_cf_principal_received)
                        ~ annual_average_temp_lag2 + annual_W_lag2 |
                          iso3c[year] + year, # trends
                        data = analysis_ccpi) %>% 
  as.list()

m_ccpi_annual3 <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                          scaled_elec_renew, scaled_fit, scaled_quota,
                          scaled_delegates, scaled_gcg,
                          scaled_cf_total_provided, scaled_cf_principal_provided,
                          scaled_cf_total_received, scaled_cf_principal_received)
                        ~ annual_average_temp_lag3 + annual_W_lag3 |
                          iso3c[year] + year, # trends
                        data = analysis_ccpi) %>% 
  as.list()

m_ccpi_season1 <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                          scaled_elec_renew, scaled_fit, scaled_quota,
                          scaled_delegates, scaled_gcg,
                          scaled_cf_total_provided, scaled_cf_principal_provided,
                          scaled_cf_total_received, scaled_cf_principal_received)
                        ~ hottest_season_temp_lag1 + season_W_lag1 |
                          iso3c[year] + year, # trends
                        data = analysis_ccpi) %>% 
  as.list()

m_ccpi_season2 <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                          scaled_elec_renew, scaled_fit, scaled_quota,
                          scaled_delegates, scaled_gcg,
                          scaled_cf_total_provided, scaled_cf_principal_provided,
                          scaled_cf_total_received, scaled_cf_principal_received)
                        ~ hottest_season_temp_lag2 + season_W_lag2 |
                          iso3c[year] + year, # trends
                        data = analysis_ccpi) %>% 
  as.list()

m_ccpi_season3 <- feols(c(scaled_climate_laws, scaled_ecp_zeros,
                          scaled_elec_renew, scaled_fit, scaled_quota,
                          scaled_delegates, scaled_gcg,
                          scaled_cf_total_provided, scaled_cf_principal_provided,
                          scaled_cf_total_received, scaled_cf_principal_received)
                        ~ hottest_season_temp_lag3 + season_W_lag3 |
                          iso3c[year] + year, # trends
                        data = analysis_ccpi) %>% 
  as.list()

list_ccpi <- do.call(c, list(m_ccpi_annual1, m_ccpi_annual2, m_ccpi_annual3,
                             m_ccpi_season1, m_ccpi_season2, m_ccpi_season3))

results_ccpi <-  setNames(as.data.frame(matrix(NA, nrow = 66, ncol = 10)), 
                            c("outcome", "treatment", "beta", "std_error",
                              "t_stat", "p_value", "nobs", "n_countries",
                              "min_year", "max_year"))

for (i in 1:66){
  results_ccpi[i, 1] <- list_ccpi[[i]]$fml %>% gsub("\\~", "", .) %>% .[2] # outcome
  results_ccpi[i, 2] <- list_ccpi[[i]]$fml %>% gsub("\\~", "", .) %>% .[3] # treatment
  results_ccpi[i, 3] <- list_ccpi[[i]] %>% tidy() %>% filter(!grepl("_W_",term)) %>% select(2) # beta
  results_ccpi[i, 4] <- list_ccpi[[i]] %>% tidy() %>% filter(!grepl("_W_",term)) %>% select(3) # std error
  results_ccpi[i, 5] <- list_ccpi[[i]] %>% tidy() %>% filter(!grepl("_W_",term)) %>% select(4) # t-stat
  results_ccpi[i, 6] <- list_ccpi[[i]] %>% tidy() %>% filter(!grepl("_W_",term)) %>% select(5) # p-val
  results_ccpi[i, 7] <- list_ccpi[[i]]$nobs # number of observations
  results_ccpi[i, 8] <- length(unique(list_ccpi[[i]]$fixef_id$iso3c)) # number of countries
  results_ccpi[i, 9] <- min(attr(list_ccpi[[i]]$fixef_id$year, "fixef_names")) # start year of analysis
  results_ccpi[i, 10] <- max(attr(list_ccpi[[i]]$fixef_id$year, "fixef_names")) # end year of analysis
}
results_ccpi$treatment <- str_sub(results_ccpi$treatment, start = 1L, end = -17L)

# results_ccpi 


## Combine results from baseline trends and CCPI models

ccpi_trends_combine <- cbind.data.frame(results_ccpi %>% 
            select(outcome, treatment, ccpi_t_stat = t_stat),
          results_trends %>% 
            select(baseline_t_stat = t_stat))

ccpi_trends_combine %>% 
  ggplot(., aes(baseline_t_stat, ccpi_t_stat)) +
  coord_cartesian(xlim = c(-5,5), ylim = c(-5,5)) +
  geom_abline(slope = 1, intercept = 0, linetype = 2, color = "grey50") +
  geom_vline(xintercept = c(-2.84, -1.96, 1.96, 2.84),
             color = c("#79b321", "#456613", "#456613", "#79b321"),
             linetype = 2,
             size = c(1.5, 1, 1, 1.5)) +
  geom_hline(yintercept = c(-2.84, -1.96, 1.96, 2.84),
             color = c("#79b321", "#456613", "#456613", "#79b321"),
             linetype = 2,
             size = c(1.5, 1, 1, 1.5)) +
  geom_point(size = 2) +
  scale_x_continuous(breaks = c(-5, -2.84, -1.96, -1, 0, 1, 1.96, 2.84, 5)) +
  scale_y_continuous(breaks = c(-5, -2.84, -1.96, -1, 0, 1, 1.96, 2.84, 5)) +
  labs(x="T-statistics in baseline model", 
       y="T-statistics in CCPI sample") +
  theme_tufte(base_family = "Helvetica") +
  theme(legend.position = "none",
        panel.background = element_rect(fill=NA),
        # panel.grid.major.x = element_line(colour = "grey90"),
        # panel.grid.minor.x = element_line(colour = "grey90"),
        # panel.grid.major.y = element_line(colour = "grey90"),
        panel.ontop = FALSE,
        text=element_text(size=16))

ggsave("text/robustness_ccpi.pdf", units = "in", height = 8, width = 10)


